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ABSTRACT 

We investigate the large-scale clustering of radio sources in the FIRST 1.4-GHz 
survey by analysing the distribution function {counts in cells). We select a reliable sam- 
ple from the the FIRST catalogue, paying particular attention to the problem of how 
to define single radio sources from the multiple components listed. We also consider 
the incompleteness of the catalogue. We estimate the angular two-point correlation 
function w(0), the variance ^2, and skewness "J 3 of the distribution for the various 
sub-samples chosen on different criteria. Both w{9) and \&2 show power-law behaviour 
with an amplitude corresponding a spatial correlation length of ro ~ 10/i -1 Mpc. We 
detect significant skewness in the distribution, the first such detection in radio sur- 
veys. This skewness is found to be related to the variance through ^3 = 83(^2)°', 
with a = 1.9 ± 0.1, consistent with the non-linear gravitational growth of pertur- 
bations from primordial Gaussian initial conditions. We show that the amplitude of 
variance and skewness are consistent with realistic models of galaxy clustering. 

Key words: galaxies: clustering - radio galaxies - large-scale structure 



1 INTRODUCTION 

Surveys of optical and infrared galaxies have revealed the 
rich structure of the universe and have provided much infor- 
mation on the large-scale structure out to redshifts z ~ 0.1 
(e.g. APM, Maddox et al., 1990 and IRAS, Fisher et al., 
1993). In contrast radio sources, although representing only 
a small fraction of all galaxies, can be detected over signif- 
icant cosmological distances (up to z ~ 4), sampling much 
larger volumes of space and therefore with the potential to 
provide information on much larger physical scales. 

It had been suggested (Kaiser, 1984) that galaxies form 
preferentially in high-density peaks of the underlying mass 
distribution. If it is true, then the statistics of galaxy distri- 
butions provide us with information, although biased, about 
the underlying matter distribution. 

It has long been debated whether radio galaxies are clus- 
tered or isotropic on the largest scale. The study by Webster 
(1976), which looked at 8000 radio sources, found < 3% vari- 
ability in the number of sources in randomly-placed 1-Gpc 
cubes. This led to a widely-accepted view that radio sources 
were isotropically distributed. Even if this were not true, the 
large range in intrinsic luminosity of radio sources might ef- 
fectively wash out structures when the distribution is pro- 
jected onto the sky with information about radial distribu- 
tion effectively lost. Other studies (Seldner & Peebels (1981), 
Shaver & Pierre (1989)) reported a detection of slight clus- 
tering of nearby radio sources, while Benn & Wall (1995) and 



Baleisis et al. (1998) discussed other measures of anisotropy 
in radio surveys. Clustering in the 4.85 Ghz Green Bank 
(87GB: Gregory & Condon, 1991) and in the Parkes-MIT- 
NRAO (PMN: Griffith & Wright, 1993) catalogues was stud- 
ied by correlation- function analysis (Kooiman et al., 1995; 
Sicotte, 1995; Loan, Wall & Lahav, 1997). These studies 
indicated that radio objects are actually more strongly clus- 
tered than local optically-selected galaxies. This conclusion 
was confirmed by correlation analysis of the FIRST survey 
(Cress et al., 1996). 

One of the possible ways of investigating clustering 
properties of radio sources is by means of the distribution 
function (counts in cells) i.e. the probability of finding N 
galaxies in a cell of particular size and shape. This analy- 
sis includes all the moments of the distribution function and 
therefore provides a more complete description of large-scale 
structure. Furthermore it can be shown (Peebles, 1980; Frie- 
man & Gaztanaga, 1994) that the higher-order moments of 
the galaxy distribution can be used as a test of non-linear 
models for large-scale structure. 

This paper presents a counts-in-cells analysis carried 
out for the FIRST radio survey. We focus on the second and 
third moment of the distribution together with the angular 
two-point correlation function of the sample, and we test 
the predictions of different cosmological models. We report 
a detection of skewness S3 = ^3/^2 = const (with <&2 and 
^3 defined as the second and the third moment of the dis- 
tribution). Our measurement accords with the hypothesis 
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of non-linear growth of observed structures by gravitational 
clustering from initially- Gaussian density fluctuations. 

In Section 2 we describe the catalogue in its original 
form and explain the procedures providing us with modified 
samples analysed in the rest of the paper. In Section 3 we 
present the results of our analysis for the angular two-point 
correlation function, and the angular second and third mo- 
ments of the distribution. Section 4 discusses the deprojec- 
tion of our 2-d measurements to estimate the quantities de- 
scribing spatial distribution. Section 5 summarises our con- 
clusions. 



artificially increase the observed clustering in the cata- 
logue. Probable side-lobes have been identified using an 
oblique decision-tree program and are flagged in the cata- 
logue (White et al., 1997). We simply rejected all the flagged 
sources from the original catalogue. We investigate the prob- 
lem of multi-component sources in some detail in the next 
section. 

We initially chose the flux limit of 1 mjy, correspond- 
ing to a 5-rms source detection (as discussed in more detail 
in Section 2.3). In this way we chose a sample of 189,689 
sources from the raw catalogue. We label this sample as CI. 



2 THE DATA 

2.1 The Public Catalogue 

The FIRST (Faint Images of the Radio Sky at Twenty cen- 
timetres) survey (Becker et al., 1995) began in the spring 
of 1993 and will eventually cover 10,000 square degrees of 
the sky in the north Galactic cap. The VLA is being used 
in B-configuration to take 3-min snapshots of 23.5-arcmin 
fields arranged on a hexagonal grid. The observations are 
at 1.4 GHz and the beam-size is 5.4 arcsecs, with an rms 
sensitivity of typically 0.14 mjy/beam. A map is produced 
for each field and sources are detected using an elliptical 
Gaussian fitting procedure (White et al., 1997). The 5-rms 
source detection limit is roughly 1 mjy. This survey is 50 
times more sensitive than any previous large-area radio sur- 
vey (see White et al., 1997 for details), leading to a high sur- 
face density of objects in the catalogue (~ 100 per square 
degree). It therefore provides an excellent tool for investi- 
gating the clustering properties of faint sources. 

The catalogue is still in the process of construction, but 
nevertheless is publically accessible. We used the 27 Feb 1997 
version which contains approximately 236,000 entries and 
is derived from the 1993 through 1996 observations cover- 
ing about 2575 square degrees, including most of the area 
7 h 20 m < RA(2000) < 17 /l 20 m , 22.2° < Dec < 42.5°. The 
sky coverage is available in the form of a map of rms sen- 
sitivity at 3-arcmin resolution. To simplify our clustering 
analysis we restricted the area to 7 h 44 m < RA < l7 h 2G m , 
22.4° < Dec < 41.8°, which has essentially complete cover- 
age. 

Within this area there are regions of low source-density 
because the original catalogue includes only sources brighter 
than 5 rms. In calculating the counts-in-cells statistics in the 
following sections we used the coverage map of the survey to 
'mask out' all those cells in which the noise was large enough 
to reduce the number of images in the catalogue, ie 0.2 mjy 
for a flux limit of 1 mjy, and 0.6 mjy for a flux limit of 3 mjy. 
When we analyse cells significantly larger than the coverage 
pixels, we reject a cell only if the fraction of useful area is 
less than 0.75. We keep cells with a ratio " ae /"' area > 0.75, 

1 total area ' 

and correct the count in each cell by the ratio of useful to 
total area. For the angular correlation-function analysis we 
corrected for these effects by generating random catalogues 
modulated by the rms sensitivity from the coverage map. 

The public catalogue is not immediately suitable for 
clustering analysis, as there are spurious images from side- 
lobes of bright sources, and single sources are frequently 
listed as multiple components. These "extra" sources would 



2.2 Multi-component Sources 

Radio sources often have widely-separated components cor- 
responding to the nucleus with hot-spots along and at the 
end of the jets. If the shape of a source is complex, the 
source detection algorithm fits several components to repro- 
duce the shape. These two effects mean that a single radio 
source can appear in the catalogue split into two or more 
sub-components. If analysed as single sources, these compo- 
nents would seriously distort clustering measurements, with 
effects particularly serious for higher order moments which 
are more sensitive to tight groupings. 

We have investigated techniques to identify tight groups 
of sources likely to be sub-components of a single radio 
source. For each group we replace all of the sub-components 
by a single source with flux equal to the sum of the individ- 
ual fluxes, placed at the mean position. As described in the 
following sections, this recombining procedure is of crucial 
importance in the analysis of the clustering. 

Following Cress et al. (1996), using a percolation tech- 
nique we first simply identified all groups of sources within 
02° of each other; each such group was replaced by a single 
source. This procedure found 25447 groups with a mean of 
2.36 sources per group. A histogram of the number of groups 
as a function of the number of sources in the group is shown 
in Figure Since the position and flux of the new compos- 
ite source is different from the individual sub-components, 
it is possible that the new source can be grouped with other 
sources. To ensure that all groups were located, we re-ran the 
group-finding procedure on the revised catalogue repeatedly 
until no new groups were found. This produced just 10 extra 
groups, leaving a catalogue of 155,084 sources. We refer to 
this sample as C2. 

This procedure efficiently combines sub-components 
into a single composite source, but can also combine 
physically-distinct sources into a single source, with dras- 
tic effect on the apparent clustering in the catalogue (see 
Figure ^). We have attempted to improve on this simple 
approach by varying the link-length in the percolation pro- 
cedure according to the flux of each source. This is motivated 
by the well-known 6 — S relation, extended by Oort (1987) 
to low flux densities, and shown to follow 6 oc yS. 

In order to determine a reasonable flux-separation re- 
lation, we plot the sum of the fluxes of the components of 
each apparent double versus their separation up to 0.05°, as 
shown in Figure ^. Superimposed on this scatter plot, we 
show contours of the density of points in the S-flux plane. 
There are two distinct regions of high density: one on the 
lower-right and the other on the lower-left of the figure. We 
inspected maps from the FIRST survey showing 4.5-arcmin 
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5 10 
number of members in group 

Figure 1. Number of component groups derived after the com- 
bining procedure vs the number of objects in each group. The 
solid line shows the results obtained by combining sources ac- 
cording to their separation (Qu n u = 0.02°); the dotted line is 
given by adding a procedure for removing spurious aggregations 
according to a flux-separation relationship, while the dashed line 
includes a constraint on the flux ratio. 




Figure 2. Sum of the fluxes of the components of double sources 
as a function of their separation. We chose to combine pairs which 
lie to the left of the dashed line (see text for details). 



squares around a sample of these objects, and found that the 
peak on the left consists mostly of sub-components of sources 
as evidenced by trailing substructures between them. The 
peak on the right consists of independent pairs of sources. 
We set the maximum link-length to be proportional to 



E 
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f I ux 1 /mJy 

Figure 3. Distribution of flux densities for members of 'double 
sources' delineated by the combining procedure described in the 
text. 



the square-root of the summed flux, Ftot, 



dlink = 100 



( Ptot Y 
\ 100 J 



(1) 



This line is shown by the dashed line in Figure |2j and roughly 
follows the minimum source-density between the two peaks. 
Varying the link-length with flux in this way combines bright 
sub-components even at relatively large separation, whilst 
keeping faint sources as single objects. As is shown in Fig- 
ure |l| (dotted line), this procedure does not seriously affect 
the number of doubles and triples, while strongly decreasing 
the number of spurious aggregations. 

As a further step in identifying groups of sources that 
should be combined, for each 'double source' identified, we 
plotted the flux of the components against each other (Fig- 
ure ^). The two fluxes are highly correlated, with most pairs 
lying within a narrow band about the value fluxl/flux2 = 
1. Such a correlation is expected for real double radio 
sources, for which the flux from the two lobes is correlated. 
Source pairs which our linking criterion does not combine 
show no correlation between the fluxes of the two compo- 
nents. This flux correlation suggests a further criterion to 
restrict the combination of sources to only physically as- 
sociated objects: we combine pairs of sources only if their 
fluxes differ by a factor less than 4 (i.e. 1/4 < < 4). 

Although somehow ad hoc, this upper limit on the flux ratio 
seems to be a reasonable guess given that it delimits a region 
on the fluxl — flux2 plane that, even though quite narrow 
(as it has to be for physically associated objects), encloses 
roughly 95% of the double sources as found through equa- 
tion jjj. Stability of clustering measurements as a function of 
the collapsing procedure will be explored elsewhere. 

The final number of groups obtained by including this 
further constraint is 17554 and their distribution as a func- 
tion of number of members is shown by the dashed line in 
Figure |l|. While keeping almost all the doubles and triples, 
the number of spurious groups formed by multicomponents 



© 0000 RAS, MNRAS 000, 000-000 



4 M. Magliocchetti, S.J.Maddox, 0. Lahav, J. V. Wall 




Log(S) 

Figure 4. Differential number counts versus the integrated flux 
density for the three catalogues CI (filled triangles), C2 (filled 
squares) and C3 (open circles), obtained from the original FIRST 
survey as explained in Sections 2.1 and 2.2. The dashed line shows 
the distribution as determined by Windhorst et al. (1985). 
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Figure 5. The angular correlation functions for CI (open 
circles), C2 (open squares) and C3 (filled circles). The error bars 
show Poisson estimates for the C3 points. The lines show w mea- 
surements from Cress et al. (1996): the solid line is for the raw 
catalogue and the dotted line for their 'collapsed' (component- 
combined) catalogue. 



is now reduced to a very small fraction (< 1CF 2 ) of the 
whole catalogue. In this way we ended up with a catalogue 
C3 formed by 167,433 sources. 

2.3 Incompleteness 

Another important issue, noted by White et al. (1997), is 
the incompleteness of the survey. Becker, White & Helfand 
(1995) estimate the catalogue to be 95% complete at 2 mjy 
and 80% complete at 1 mjy. This could affect our results 
in two ways. First, if there are any spatial inhomogeneities 
in completeness, spurious clustering will be introduced. Sec- 
ond, even if the incompleteness is uniform, the redshift dis- 
tribution will be different from a complete flux-limited sam- 
ple, and different in an unpredictable way. This difference 
will change the relation between 2-d and 3-d quantities, and 
so bias the 3-d measurements, as discussed in Section 4. 

In Figure ^ we plot the differential source counts nor- 
malised to the distribution expected for a non-evolving pop- 
ulation in Euclidean space vs the integrated flux density for 
the three catalogues CI (filled triangles), C2 (filled squares) 
and C3 (open circles), obtained from the original FIRST 
survey as described in Sections 2.1 and 2.2. The dashed line 
shows the power-law fit to the 4£L distribution determined 
by Windhorst et al. (1985). There is a drop in the number 
counts for fluxes smaller that 3 mjy, and this lack of faint ob- 
jects becomes very important below 2 mjy. Furthermore the 
summation of fluxes in the recombination of sub-components 
affects the incompleteness of the resulting catalogues. This 
is especially important in the faint-flux region. Note that as 
we chose an initial threshold of 1 mjy, we may lose pairs 
with f\ + /2 > 3 mjy, but with either fx or fa < 1 mjy. 

From Figure ^, it is a reasonable assumption that the 
survey is complete at fluxes greater than 3 mjy. Hence we 



set a flux limit of 3 mjy for the catalogues used in our anal- 
yses. At this flux limit the C1,C2 and C3 samples contain 
respectively 101787, 83287 and 86074 objects. In the rest of 
the paper we adopt the C3 version of the catalogue, but we 
compare to the results with those obtained for CI and C2. 



3 CLUSTERING ANALYSIS 

3.1 The Angular Two- Point Correlation Function 

Correlation-function analysis has become the standard way 
to quantify the clustering of different populations of astro- 
nomical sources. Specifically the angular two-point correla- 
tion function Wx2 = w(8) gives the excess probability, with 
respect to a random Poisson distribution, of finding two 
sources in the solid angles SQx Sfl2 separated by an angle 6, 
and it is defined as 

5P = nSVLxS^ [1 + w(6)] (2) 

where n is the mean number density of objects in the cata- 
logue under consideration. 

We measured w for the three catalogues CI, C2 and C3 
using various techniques. First we generated a random cata- 
logue of positions and fluxes and selected sources above the 
sensitivity limit from the coverage map. Then we counted 
the number of distinct data-data pairs (DD), data-random 
pairs (DR), and random-random pairs (RR) as a function 
of angular separation. We then measured w from the three 
commonly-used estimators (Hamilton, 1993; Landy & Sza- 
lay, 1993; Peebles, 1980), 



DD+RR-DR. 
~ DR 
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2DD 
DR 



(3) 



We also counted the number of sources in equal-area cells 
as discussed in section 3.2 and estimated w from the cell 
counts, 



(N1N2) 
(N 1 }(N 2 ) 



- 1 



(4) 



For all catalogues, wis and Wh were essentially identical, 
and not significantly different from w ce u. The estimate Wdr 
agrees perfectly on scales 9 < 1°, but is slightly lower at 

In Figure ^ we show Wh for the three catalogues CI, 
C2 and C3; the error bars show Poisson estimates for the 
C3 points. As the distribution is clustered these estimates 
only provide a lower limit to the real errors (see e.g. Mo et 
al., 1992) that might be obtained for instance by using the 
bootstrap re-sampling method. More precise estimates will 
be given in section 3.2 when we will move on to more quan- 
titative results. 

Note that, as it will be extensively discussed later in the 
paper, versions of the same catalogue obtained by adopting 
different collapsing procedures lead to different estimates of 
the angular correlation function; as it can be seen by com- 
paring the C2 measurement with those given by the analysis 
of the CI and C3 catalogues, it turns out that the higher 
is the number of objects that have been combined together, 
the shallower is the slope of the corresponding Wh- In more 
detail the CI measurement shows a steepening at 9 < 0.02° 
caused by the over-counting of multiple-component sources. 
At very small angles (9 < 0.002°) the amplitude drops to —1, 
as close pairs of sources become unresolved by the survey. 
The C2 measurement is —1 for 9 < 0.02° because all DD 
pairs have been removed from the catalogue by the recom- 
bining procedure. Our C2 measurement should be compara- 
ble to the w estimate from Cress et al. (1996). In fact our 
measurements for both CI and C2 are in very good agree- 
ment with Cress et al. (1996) for 9 < 0.1°, but at larger 
scales we find a slightly steeper slope, leading to a difference 
of a factor of 2 by 9 ~ 2°. The discrepancy is marginally 
significant given the estimated errors, and is probably due 
to using the more recent version of the catalogue. 

The C3 measurement continues approximately as a 
power law to the limit of the survey resolution, 9 ~ 0.002°. 
This gives us some confidence that our multi-component re- 
combining procedure is reliable. 



3.2 Variance 

The galaxy distribution function (counts in cells) gives the 
probability of finding N galaxies in a cell of particular size 
and shape. In principle the counts-in-cells distribution is 
straightforward to compute; one divides the space into cells 
of a given volume (or in our case, of given area) and shape, 
and then counts the number of galaxies in each cell. Obvi- 
ously there cannot be less than zero galaxies in a cell, but 
the number of objects in each cell can grow with no upper 
bound (unless some negative feedback process takes place). 
If we define the kih moment of the counts as 

^ k = {(N-N) k ) (5) 
where N = nfl is the mean count in the solid angle Q, it is 
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Figure 6. The normalised variance a 2 vs the cell size © for our 
C3 catalogue of the FIRST survey. Errors are estimated from the 
variance in four random subsets. 



then possible to write these moments in terms of the n-point 
correlation functions (see e.g. Peebles, 1980). The second 
moment of the galaxy distribution function is related to the 
two-point correlation function W12 through the expression 



fX2 = N + (JVf*2 



(6) 



where N is the shot noise resulting from the discrete nature 
of the sources (Poisson noise), and 



* 2 



n 2 



W12 dfli dfl2 



(7) 



is the normalised variance in terms of the two-point correla- 
tion function integrated over a cell of area f2 and particular 
shape. 

By assuming a power-law form 
w{9) ^Ad 1 ' 1 , 



(8) 



and by considering square cells of size 11 = 6x6 square de- 
grees, we can evaluate the integral in equation Q (see Totsuji 
& Kihara, 1969), obtaining 



^2 



N 



(N /CI) 2 



At 



y dn 1( in2 



A 5 — I 



0) 



where C 7 is a coefficient depending on 7 which can be eval- 
uated numerically by Monte Carlo methods (see e.g. Lahav 
& Saslaw, 1992). It is therefore possible to use the a 2 - 6 re- 
lation to evaluate the two parameters (^4, 7) which describe 
the correlation function. 

To carry out our analysis we first project the coordi- 
nates of the objects in the survey onto an equal-area projec- 
tion. We then divide the region into square cells of constant 
x and y in rectangular coordinates for the projection, and 
count the number of objects in each cell. We use the cover- 
age map to discard all cells which are partially occupied, as 
described in section 2.2. The procedure is repeated for dif- 
ferent cell sizes, from Ax = Ay = O = 0.3° up to 3°; at each 
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Figure 7. The normalised variance a 2 vs cell size © for the CI 
catalogue of the FIRST survey. Errors are estimated from the 
variance in four random subsets. 
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Figure 8. The normalised variance a 2 vs cell size © for the C2 
catalogue of the FIRST survey. Errors are estimated from the 
variance in four random subsets. 



cell size, the variance a 2 is calculated from the estimator in 
equation [| 

Figure [i| shows a 2 as a function of O for the C3 cata- 
logue. The slope and the intercept of the plot are estimated 
by a least-squares procedure minimising the quantity 



y-N I log(<r)» - a- &log(8)j 

i=l 



(10) 



with b = (5 — 7), a = log(ylC 7 ). The errors Aj are obtained 
using the 'partition bootstrap method' in which the nor- 
malised variance is calculated for four subdivisions of the 
survey region and the standard deviation of these measure- 
ments at each angle is used as a measure of the error. Note 
that this is not strictly a x 2 statistic, since the individual 
points are not independent. From this analysis we find 



7 = 2.50 ±0.1 ; A = (1.06 ± 0.1) • 10" 



(11) 



where C 1 = 6.52 has been obtained by solving the integral 
in equation |^ with the condition 6 > 0.01°, appropriate to 
the component-combining algorithm used to construct C3. 

The amplitude A is smaller by about one order of mag- 
nitude than the values obtained from optical surveys (e.g. 
A = (3.82 ±0.12) ■ 10" 2 for the APM survey (Maddox et al., 
1990)); this is due to the 'washing out' of structure by the 
wide redshift span of radio surveys. The value we find for 7, 
high in comparison with optical surveys (e.g. 7 ~ 1.8, Pee- 
bles (1980), Maddox et al. (1990)), is discussed in Section 
4. 

By repeating the analysis for the CI and C2 samples 
(Figures |^ and |8), we find respectively: 

7 = 2.68 ±0.07 ' 1 - '1 " -L n rurt m-3 



A = (1.52 ± 0.06) • 10" 



and 



7 = 2.23 ±0.1 ; A = (1.42 ± 0.1) ■ 10" 



(12) 



(13) 



from equation |9j by imposing the condition 9 > 0.02°, again 
from the nature of the component-combining criterion. The 
values in equation [H] are in close agreement with those ob- 
tained by Cress et al. (1996). 

Comparison of the slope 7 in equations 0, and P 
shows how its value decreases according to the number of 
source-components combined in the different catalogues (7 
largest for no combining) . The presence of more 'sources' or 
'source components' in close proximity, regardless of their 
physical association, results in a stronger correlation on 
smaller scales. 



3.3 Skewness 

Further information on clustering is contained in the higher 
moments of the distribution, such as the skewness and the 
kurtosis. Assuming Gaussian primordial perturbations and 
linear theory, both the skewness and the kurtosis remain 
zero (Peebles, 1980). A detection of a non-zero value for 
these quantities may then be evidence of non-linear gravi- 
tational clustering. Departures from Gaussian distributions 
might also reflect non-Gaussian initial conditions predicted 
in some cosmogonies such as cosmic-string and texture mod- 
els. Detection of non-zero values for skewness and kurtosis 
is therefore of fundamental importance. 

Kurtosis is not considered here because its estimated 
value is dominated by noise. Skewness is actually observed 
in each of the catalogues CI, C2 and C3 under consider- 
ation, starting from angles Q ~ 1°. As seen in Figure 
the galaxy distribution functions are noticeably skewed to- 
wards larger values of source-numbers, the tail becoming 
increasingly prominent as the cell-size becomes smaller. To 
illustrate, Figure [lo] shows the comparison between sam- 
ple C3 and a Poisson distribution. K-S tests show that the 
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Figure 9. Distribution function for the FIRST survey for differ- 
ent cell sizes. The solid lines show the results obtained for CI, 
and the dotted and dashed lines are for C2 and C3 respectively. 
The x axes have been rescaled according to x' = x — x m , where 
Xm is the median. 
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Figure 10. Distribution function for the FIRST survey for dif- 
ferent cell sizes. The solid line shows the results for sample C3 
compared to a Poisson distribution (dotted line). 



probabilities that the histograms are consistent with Poisson 
distributions are < 0.01%. 

Using an approach similar to that described in Section 
3.2, the third moment (skewness) is related to the three- 
point correlation function u)\23 through the expression 



^ 3 = JV + 3(jV) 2 # 2 + (.W) #3, 



(14) 



-2.5 - 



1.5 



A 



_L 



_L 



1.4 



-1.6 -1.8 -2 
log(* a ) 

Figure 11. The normalised variance ^2 = £2 2 <r 2 vs the nor- 
malised, dimensionless skewness \E"3 for C3. Errors have been cal- 
culated as in Section 3.2. 



where 



#3 



1 



W123 dQi dQ,2 dQ,3 



(15) 



is the normalised skewness in terms of the three-point cor- 
relation function integrated as in equation [?]. 
An estimator for the normalised variance is 



and an estimator for the normalised skewness is 
7V-3(iV) 2 *2 



*3 = 



fJ-3 



(JV) 8 



(16) 



(17) 



It was shown (Peebles, 1980; Juszkiewicz et al., 1995; Coles 
& Frenk, 1991) that if the initial perturbations were Gaus- 
sian and they grew with cosmic time purely due to gravity, 
then in second-order perturbation theory (i.e. in the quasi- 
linear regime, 8 = — ~ 1) 



SS(R) 



(S 3 )r 



f~(n + 3) 



(18) 



where S3 (R) is the spatial normalized skewness in randomly 
placed spheres of radius R and n is the index of the primor- 
dial power-spectrum of fluctuations. It follows that in this 
case the projected quantities obey ^3 oc ($2)", with a — 2. 
For local surveys (e.g. Peebles, 1980), a — 2 is also expected 
from the empirical relation between the 3-point and 2-point 
correlation functions: 



C(ri,r a )=0K(ri)f(r 3 )+«r 1 )f(ri2) + 
+£(r 2 )£(ri 2 )](xe 2 (r), 



(19) 



where £(n,r2) is the spatial three-point correlation func- 
tion. 

Figure |ll] shows log^s) vs log($2) for angles 0.3° < 
< 1°; a least-squares fit relates the quantities by the form 



*3 = S 3 (0) *2 



(20) 
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Figure 12. S3 as a function of for C3, assuming a = 2. Errors 
are estimated as in Section 3.2. 
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Figure 14. S3 as a function of © for CI assuming a = 2. Errors 
are estimated as in Section 3.2. 
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Figure 13. The normalised variance ^2 = fl 2 A vs the nor- 
malised, dimensionless skewness "1/ 3 for CI. Errors have been cal- 
culated as in section (3.2). 
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Figure 15. The normalised variance ^2 = f2 2 cr 2 vs the nor- 
malised, dimensionless skewness ^3 for C2. Errors have been cal- 
culated as in Section 3.2. 



with a — 1.9±0.1 and Sj^S) = 9±3 and errors estimated as 
in Section 3.2. Therefore the value found for a supports the 
assumption of gravitational growth of perturbations from 
Gaussian initial conditions. 

To test the constancy of S3 as predicted by the hierar- 
chical theory, we plot this quantity as a function of the cell 
size on the assumption of a — 2 (Figure ^). The plot 
supports the hypothesis of constant Ss(0) = 12.5 ± 1.5. 

The same analysis for the catalogue CI (containing all 
sources), and C2 (generated by combining 'sources' only ac- 
cording to their separation), finds: a = 2.0 ± 0.1, Ss(Q) = 
4 ±2; S3 (9) = 4.7 ±0.3 for CI (see Figures [[3] and [jj) , and 



a = 1±2, S3 (6) = 1±3; S 3 (6) = 17±5 for C2. In this latter 
case the estimate of the quantities is dominated by noise (see 
Figures [l5] and [l6|). The results are highly sensitive to the 
process of 'combining' multi-component sources. Combin- 
ing two or more 'sources' in close proximity strongly affects 
the distribution function in two ways. It removes the tail, 
thereby reducing skewness, and at the same time it reduces 
the variance as more cells will contain numbers close to the 
average of the distribution. The effects become rapidly more 
important as the number of combined sources increases; this 
explains the remarkable difference in the values obtained 
from the three catalogues. In particular the results for the 
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Table 1. COUNTS IN CELLS ANALYSIS RESULTS 





Cl 


C2 


C3 


7 


2.68 ±0.07 


2.2 ±0.1 


2.50 ±0.1 


A 


(1.52 ±0.06) • 10" 3 


(1.4 ±0.1) ■ 10~ 3 


(1.06 ±0.1) ■ 10" 3 


a 


2.0 ±0.1 


1±2 


1.9 ±0.1 


Sa(6) 


4±2 


1±3 


9± 3 


5 3 (e) Q=2 


4.7 ±0.3 


17 ±5 


12.5 ± 1.5 
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Figure 16. S3 as a function of © for C2 assuming a = 2. Errors 
are estimated as in Section 3.2. 



C2 catalogue show that the variance and the skewness of its 
distribution become so small as to be dominated by noise. 
Note that the differences between C2 and C3 are due to less 
than 3000 'objects' (~ 4%). 

Table 1 summarises the results obtained from the 
counts-in-cells analysis for the three catalogues Cl, C2 and 
C3. 



4 RELATION TO SPATIAL QUANTITIES 

4.1 The Spatial Correlation Function 

Once both the cosmological model and the selection function 
of the sample are specified, the standard way of relating the 
spatial (3d) correlation function £(r,z) to the angular (2d) 
one w(9) is via the relativistic Limber equation (Peebles, 
1980). For Einstein-de Sitter universe (fi = 1,A = 0), 



w{9) 



ir JT z 4$2 (aOf (r, z)dx du 



where x is the comoving coordinate 



(21) 



(22) 



and the selection function &(x) satisfies the relation 



N - 



&(x)x 2 dx ■■ 



1 



N(z)dz, 



(23) 



10 "" s Jo 

in which M is the mean surface density on a surface of solid 
angle £l 3 and N(z) is the number of objects in the given 
survey within the shell (z, z + dz). 

If we assume for £ (r, z) a power-law redshift- dependent 
form 



(1 + *) 



-(3+£) 



(24) 



7 constant with z, where r is the proper coordinate, ro the 
correlation scale length at redshift z = and e a parameter 
describing the redshift evolution of the spatial correlation 
function, then, in comoving coordinates r c = r(l + z), £ 
assumes the form 



(1 + 



\7-(3+e) 



(25) 



Specific values for e have the following significance: e = 
implies constant clustering in proper coordinates; e = 7 — 3 
implies constant clustering in comoving coordinates; and e = 
7 — 1 represents growth of clustering under linear theory 
(Peebles, 1980; Treyer & Lahav, 1996). 



In the small-angle approximation 



(u a + 



x 2 9 2 ) 1 ^ 2 /(I + z) we then find that the amplitude A defined 
in equations ^ and|^ is given by the expression (e.g. Loan et 
al., 1997; Peebles, 1980): 



r(l)r(^i) 



(26) 



~ e dz 



j-[l-(l + z)-^N 2 {z)(l + z) 
[/" N(z)dz] 2 

where P is the Gamma function. This expression, consider- 
ing the functional form (^) for w(8), directly fixes the value 
of ro once the evolution parameter e is given. 



4.2 The Spatial Skewness 

An analysis analogous to that in Section 4.1 provides the 
3-dimensional value for the skewness S^(R) in randomly- 
placed spheres of comoving radius R. This is related to the 
projected (2-dimensional) quantity §3(6), in the case of cir- 
cular cells of radius 0, through the equation (Gaztahaga, 
1995) 



s 3 (e) = r 3 s* 3 (R)(^) 



(27) 
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Redshift 

Figure 17. Redshift distribution N(z) for the radio source pop- 
ulation at 1.4 GHz at a flux limit of 3 mjy. The dotted curves 
represent the 6 models (1-4, 6-7) taken from Dunlop & Peacock 
(1990); the solid curve is the average. 



Redshift 

Figure 18. Redshift distribution N(z) for the radio source pop- 
ulation at 1.4 GHz at a flux limit of 3 mjy. The dotted lines are 
models 6 and 7 from Dunlop &; Peacock (1990); the solid line is 
the average. 



where Cz ~ 1 and B3 ~ 1 (see Gaztanaga 1994 for further 
details), and r$ can be expressed as 



with 

/*oo 

I k = x^x^x^-^-^il + zf^-^ 1 -^, (29) 
Jo 

<&(ie) being defined by equation 

4.3 The redshift distribution N(z) 

To obtain both the spatial correlation function and the spa- 
tial skewness from the respective projected quantities, we 
need N(z), the redshift distribution for radio sources. There 
is not even a completely-identified sample of radio sources 
at 3 mjy, leave alone one with measured redshifts. How- 
ever, the Dunlop & Peacock (DP, 1990) models of epoch- 
dependent luminosity functions for radio sources can be used 
to make estimates of N(z). These models were derived with 
Maximum Entropy analysis to determine the coefficients 
for polynomial expansions to represent the epoch-dependent 
luminosity function; the approach incorporated the then- 
available identification and redshift data for complete sam- 
ples from radio surveys at several frequencies. From this we 
have adopted three different models of N(z) for subsequent 
deprojection analysis. 

The first model (Ni (2)) is shown in Figure |l7]. The dot- 
ted lines represent the six models (1-4, 6-7) taken from DP 
and the solid line, the average, is our Ni(z). The large spread 
indicates the uncertainty due to incomplete or statistically- 
limited redshift data as well as the extrapolation of the data 
to such a low flux density. Nevertheless there are two in- 
teresting features which are general. The spike seen at small 



redshifts indicates that at such low flux densities, the lowest- 
power tail of the local radio luminosity function begins to 
contribute substantial numbers of low-redhift sources. This 
spike is almost certainly underestimated; the DP analysis 
did not encompass the evolving starburst-galaxy population 
(e.g. Windhorst et al., 1985) now believed to constitute a 
majority of sources at mjy levels (see Wall, 1994 for an 
overview). The second feature is the prominent maximum 
around z ~ 1 displayed by all models, confirming that the 
median redshift for radio sources in radio surveys over a 
wide flux-density range is a factor of 10 higher than that of 
wide-field optical surveys. 

The second model (N2(z)) considered for our deprojec- 
tions is represented by the solid line in Figure |l^ and is 
obtained by averaging models 6 and 7 from DP, models of 
pure luminosity evolution. The more prominent spike at low 
redshifts may be a more realistic representation of the total 
population. 

The third model (Ns(z)) is simply from Ni(z) with the 
small-z spike patched out. 

We take Ni(z) as our best estimate, and use ^2(2) and 
Ns(z) to test how sensitive the deprojection is to the form 
of N(z). In particular we wish to see how the results for the 
correlation length ro and the quantity r% in equation ^ are 
affected by "extreme" models (one dominated by the spike 
and one completely without it). 

4.4 The Predicted Angular Correlation Function 

The values obtained for the slope 7 in Section 3.2 are 
high compared to the canonical value of 1.7 measured for 
optically-selected galaxies. This may be due to the intrinsic 
clustering properties of radio galaxies, but in this section we 
consider another possibility. The depth of the survey, z ~ 1, 
means that the angular scales 0.3° < 8 < 3°, where we 
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Figure 19. CDM prediction for the correlation function with 
T = 0.2. The y axis is the absolute value of For r ~ lOOh," 1 
Mpc, £ (r) becomes negative, and we plot its absolute value as the 
dotted line. 



0.01 



0.001 



0.0001 



10" 



~1 1 1 — I — I — I I I 



~i 1 r 



0.1 



e 



Figure 20. CDM predictions for w(8), showing the effects of 
varying the cosmological parameters, with N(z) = Ni(z) as de- 
scribed in Section 4.3. The dotted line is for f2 = 1, h = 1, the 
solid line for Q = 0.3, h = 0.65, and the dashed line for Q = 0.5, 
h = 0.5. 



measure clustering, correspond to linear scales large enough 
such that the correlation function £ (r, z) is no longer a simple 
power-law, but has steepened from that given by 7 ~ 1.7. 

If this is the case, choice of which value of 7 to use in the 
deprojections of equations ^ and ^is not straightforward. 
We have measured the slope 7 = 2.50 at large scales (Section 
3.2), but simply using the same slope for all scales would lead 
to an overestimate of w(9) and therefore of £(r, z) (through 
equation ^) at small scales. Both ro and are very sensitive 



to 7 (see below) and therefore the assessment of its proper 
value is very important for our conclusions. 

To investigate what values are expected for 7 we have 
taken a theoretical correlation function £(r, z) up to scales 
r ~ lOOO/i" 1 Mpc and projected it to get w(8) by directly 
integrating equation [2lJ By doing the integrations numeri- 
cally we do not make the assumption that £ is a power law, 
and do not use any of the approximations used in Section 
4.1. In Figure ^ we plot the £(r) used. On large scales it 
is the linear prediction for a CDM model with F = 0.2. 
On small scales the linear-theory prediction underestimates 
the true amplitude, and as a rough approximation to allow 
for this, we simply extrapolated a power-law of slope —1.7 
for scales r < 10/i -1 Mpc, as observed from the APM cor- 
relation function (Maddox et al., 1990). Given the simple 
heuristic aims of this section, and the large uncertainties 
in the evolution of galaxy clustering and bias, more sophis- 
ticated methods to account for non-linearity (e.g. Peacock 
& Dodds, 1996) are unnecessary, though we intend a more 
careful study in a future paper. The overall normalization is 
set so that ro = 5Ah~ 1 Mpc. 

Given this form for £(r), the shape and the amplitude 
of w(0) still depend on the functional form of N(z) and the 
values of S7, h. Figure ^ shows w{6) obtained using the red- 
shift distribution Ni(z) as introduced in Section 4.3, with 
three different sets of cosmological parameters. The curves 
are significantly different, but they all show the same power- 
law behaviour for small angles, with a slope of ~ —0.7, 
corresponding to the expected (1 — 7) with 7 = 1.7; then, 
at 6 ~ 0.3°, there is a break from this power law towards 
steeper slopes, corresponding to 7 ~ 3.5 for 9 > 2°. 

In the interval 0.3° 6 3° the curves are close to a 
power-law with slope ~ —1.5. This is in agreement with the 
value of 7 w 2.5 obtained in Section 3.2. The smaller the 
value of f2, the steeper is the effective slope. Of the three 
curves in the range 0.3° 5s 8 5s 3°, that for Q. = 0.3 and 
h = 0.65 fits the data best. However this rather superficial 
analysis assumes a fixed f(r) and N(z), both of which will 
significantly change the shape of w. We intend to make a 
more thorough comparison between models and the data in 
a future paper. 



4.5 Results 

We now estimate the correlation length ro and the spatial 

normalized skewness S3 (R) using equations 

with different redshift distributions, and a range of values for 

7 and e. We confine attention to our most reliable clustering 

estimates, namely those obtained in Section 3.2 from the C3 

catalogue. 

Figure ^ shows the trends for both the correlation 
length ro and the quantity rz, as functions of the evolu- 
tion parameter e, obtained for different values of 7 and the 
three different models for the redshift distribution N(z) in- 
troduced in Section 4.3. If we fix a value for 7 and vary 
the form of N(z), the variations for ro (Sj 10%) are not as 
dramatic as those for r-j. In general, the stronger the low- 
redshift spike in N(z) the smaller the value obtained for the 
correlation length ro. If instead we fix the functional form 
for N(z), varying 7 strongly affects both ro and r^, - Note also 
that the value of the amplitude A used in the deprojection 
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Figure 21. Correlation length ro and rz ~ S3 (©)/S| (-R) as a function of the evolution parameter e for different values of 7 and 
respectively JVi, N2, N3, in the case of the C3 catalogue; the solid line corresponds to to 7 = 1.8, the dotted line to 7 = 2.0, the 
short-dashed line to 7 = 2.2, and the long-dashed line to 7 = 2.50. 



equation |2£] to obtain the correlation length ro also depends 
on 7 through the quantity C 1 (equation H). 

Armed with these results we face choosing the value for 
7 to provide the best estimates for ro and rz- Our counts- 
in-cells analysis carried out for scales Q < 0.3° is dominated 
by the shot (Poisson) noise, so we must consider cluster- 
ing on larger scales. For 0.3° < 6 < 3° we observe a slope 
7 = 2.50 but from Section 4.4 it follows that clustering on 
these angular scales is determined by the spatial correlation 
function on scales larger than lO/i" 1 Mpc where a power 
law is not a good approximation. As mentioned in section 
4.4, using 7 = 2.50 in equation [§] will overestimate w(6), 
and hence £(r) at small scales where the spatial correlation 
function is noticeably greater than 1 (see Figure |l^). This 
will severely bias the value of ro obtained from integrat- 
ing equation ^(j. This problem is even worse when we try 
to deproject the skewness, for the quantity rz appearing in 
equation Bq depends on the square of the correlation func- 



tion (see Section 3.3). The choice of 7 is further complicated 
by the fact that N(z) suggests the presence of at least two 
populations of radio sources, which may well have different 
clustering properties. Yet another uncertainty is introduced 
by fact that 7 may well be a function of redshift, even for a 
single population. 

For convenience of comparison with other studies we 
quote our estimates using the "standard value" 7 = 1.8 and 
assume stable clustering (e = 0). Adopting Sz(Q) = 12.5, 
we then find 



r ~ 9.7k- 1 Mpc; S 3 * (i?) ~ ~ 4.2, 



Ni(z) 

N 2 {z): r ~ 9.7ft" 1 Mpc; Sz{R) 
Nz(z) 



rz 
rz 



1.2, 



r ~ Wm^Mpc; S 3 *(i?) ~ ~ 8.9. 



rs 



(30) 
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We note that the values for ro are little affected by the choice 
of N(z). Our estimate of ro for radio sources is larger than 
the value for optical (ro ~ Mpc) and IRAS galaxies 

(ro ~ 4ft,— 1 Mpc) and smaller than for the cluster-cluster 
correlation length (ro ~ 14 — 20h~ 1 Mpc). 

On the other hand, we see that the spatial normalized 
skewness Ss(R) is sensitive to N(z). (We have neglected the 
small corrections for the difference between square and cir- 
cular cells.) Apart from the observational uncertainties in 
estimating Ss(R) for the radio sources, we should remem- 
ber that comparison with the theoretical value of S3 (R) for 
the mass perturbations (equation |l8| ) depends on how radio 
galaxies are biased relative to the mass distribution. This 
bias may also be epoch-dependent (e.g. Fry, 1996). Neverthe- 
less, our estimates are in accord with the prediction (equa- 
tion [t|) of S£(R) ~ 2.9 for power-spectrum P{k) ~ fc" 1 
expected (based on other observations and models) on the 
weakly non-linear scales R ~ 5 — 15ft." 1 Mpc probed by our 
analysis. (Note that angular scale of 1 degree at the median 
redshift of the survey, z ~ 1, corresponds to comoving dis- 
tance of 30/i _1 Mpc for an Einstein - de Sitter universe.) 
It is also interesting to note that our values for S3 bracket 
the values derived from the optical APM survey (S3 ~ 3.2; 
Gaztahaga et al., 1994) and from the IRAS survey (S3 ~ 2.8; 
Kim & Strauss, 1998), even though the observed difference 
in the spatial skewness could reflect differences in cluster- 
ing properties between radio sources (possible mix of more 
objects like bright ellipticals and starbursting galaxies) and 
optical or IRAS galaxies. 



5 CONCLUSIONS 

By investigating the distribution function (counts in cells) 
for radio sources of the FIRST survey, we have shown how 
to infer the clustering properties of host radio galaxies. We 
considered three different catalogues above 3 mjy, generated 
from the original survey by following different procedures for 
combining source components. Focusing on the catalogue 
obtained with the most sensible procedure, by relating the 
second moment of the distribution to the angular two-point 
correlation function, we find a power-law behaviour for w(8). 

From the analysis of the third moment we have shown 
how variance and skewness are related through the func- 
tional form ^3 = Ss('i>2) a with a ~ 2 and S3 = const with 
angular scale. By inverting the projected quantities we have 
estimated the spatial correlation length ro ~ 10ft _1 Mpc and 
the spatial skewness S-j(R) ~ 1 — 9. While the value for ro 
is relatively independent of the functional form of the red- 
shift distribution N(z), the value for skewness is strongly 
dependent on it. 

Our results indicate that the large-scale clustering of 
radio sources is in accord with the gravitational instabil- 
ity picture for the growth of perturbations from a primor- 
dial Gaussian field. Our measurements of deviations from 
a Gaussian distribution do not seem to require initial non- 
Gaussian perturbations (e.g. cosmic strings or texture). 

There are crucial observations required to further this 
analysis. We need a better understanding of populations and 
source structures at mjy levels to remove statistical uncer- 
tainties due to multi-component sources. Even more impor- 
tant is the observational measurement of N(z), from iden- 



tifications and redshift measurements for a complete sam- 
ple. Radio morphologies and optical studies for a small flux- 
limited sample would achieve both goals. 
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